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O ■ ABSTRACT 

CN ! Smoothed particle hydrodynamics (SPH) employs an artificial viscosity to properly capture 

hydrodynamical shock waves. In its original formulation, the resulting numerical viscosity 
is large enough to suppress structure in the velocity field on scales well above the nom- 
inal resolution limit, and to damp the generation of turbulence by fluid instabilities. This 
could artificially suppress random gas motions in the intracluster medium (ICM), which are 
driven by inf ailing structures during the hierarchical structure formation process. We show 
. that this is indeed the case by analysing results obtained with an SPH formulation where 

■ an individual, time-variable viscosity is used for each particle, following a suggestion by 
i/^ \ [Morris & Monaghan ( 1997). Using test calculations involving strong shocks, we demonstrate 

■ that this scheme captures shocks as well as the original formulation of SPH, but, in regions 

away from shocks, the numerical viscosity is much smaller. In a set of nine high-resolution 
simulations of cosmological galaxy cluster formation, we find that this low-viscosity for- 
mulation of SPH produces substantially higher levels of turbulent gas motions in the ICM, 

O \ reaching a kinetic energy content in random gas motions (measured within a IMpc cube) 

of up to 5%-30% of the thermal energy content, depending on cluster mass. This has also 
^ ■ significant effects on radial gas profiles and bulk cluster properties. We find a central flat- 

' tening of the entropy profile and a reduction of the central gas density in the low-viscosity 

scheme. As a consequence, the bolometric X-ray luminosity is decreased by about a factor 
of two. However, the cluster temperature profile remains essentially unchanged. Interestingly, 
\ this tends to reduce the differences seen in SPH and adaptive mesh refinement simulations of 

^ ■ cluster formation. Finally, invoking a model for particle acceleration by MHD waves driven 

by turbulence, we find that efficient electron acceleration and thus diffuse radio emission can 
be powered in the clusters simulated with the low viscosity scheme provided that more than 
5%-10% of the turbulent energy density is associated with Fast Magnetosonic Modes. 
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1 INTRODUCTION 

In hierarchical cold dark matter cosmologies, l arge structures form 
through the accretion of smaller structures fe.g. JWhite et all 19931) . 
In particular, mergers and infall of halos play a fundamental role 
in determining the structure and dynamics of massive clusters of 
galaxies, where mergers can induce large-scale bulk motions with 
velocities of the order of ^ 1000 km s~^ or larger. This results in 
complex hydrodynamic flows where most of the kinetic energy is 
quickly dissipated to heat by shocks, but some part may in principle 
also excite long-lasting turbulent gas motions. 

Numerical simulations of merging clusters (g-g-? 
ISchindler & MuelleJl99l:lRoettiger et alll99llRicker & Sarazid 
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I2OOII : iTakizawal l2005h provide a detailed description of the 
gas-dynamics during a merger event. It has been shown that 
infalling sub-clusters can generate laminar bulk flows through the 
primary cluster and inject turbulent eddies via Kelvin-Helmholtz 
(KH) instabilities at interfaces between the bulk flows and the 
primary cluster gas. Such eddies redistribute the energy of the 
merger through the cluster volume in a few turnover times, which 
corresponds to a time interval of order 1 Gyr. The largest eddies 
decay with time into more random and turbulent velocity fields, 
eventually developing a turbulent cascade with a spectrum of 
fluctuations expected to be close to a Kolmogorov spectrum. 

Observationally, spatially resolved gas pressure maps of the 
Coma cluster obtained from a mosaic of XMM-Newton observa- 
tion have indeed revealed the signature of mildly s upersonic turbu- 
lence, at least in the central regions of the cluster iSchuecker et alJ 
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l2004 . It has also been suggested that the micro-calorimeters on- 
board of future X-ray satellites such as ASTR0-E2 should be 
able to detect the turbulent broadening of the lines of heavy ions 
in excess of the thermal broadening (Inogamov & Sunvaev 2003), 
which would represent a direct measurement of cluster turbulence. 

Cluster turbulence could in principle store an appreciable frac- 
tion of the thermal energy of massive clusters, which would make 
it an important factor for understanding the structure of the ICM. 
Shear flows associated with cluster turbulence and the resulting 
dynamo processes could c onsiderably amplify the m agnetic field 
strength in the ICM (e.g.. ^ Dolag et alJll999L 12002^. In addition, 
magnetohydrodynamic waves can be efficiently injected in the ICM 
directly by shocks, by Kelvin-Helmholtz or Rayleigh-Taylor insta- 
bilities, or by the decay of turbulent eddies at larger scales. These 
waves, as well as shocks, can efficiently accelerate supra-thermal 
particles in the ICM to higher energies. Although there is still some 
debate concerning the detailed mechanism responsible for the ori- 
gin of relativi stic particles and magnetic fields in the ICM (e.g., 
lBmnettil200il the presence of relativistic electrons and of ~ jiG- 
strength magnetic fields in the ICM is proven by non-thermal emis- 
sion studied with radio observations and possibly observations of 
hard X-Ray emission (e.g., Feretti et al. 2002; Fusco-Femiano et al. 
I2OO3L for a review). In addition, the occurrence of non-thermal phe- 
nomena is found to be re lated to th e dynamical state and mas s of the 
parent cluster ( Giovanni ni et alil999: Buote 2001: Schueck er et alJ 
l200ll: lFeretti 2002), which suggests a connection between cluster 
mergers and non-thermal activity. 

Despite this potentially significant relevance of turbulence for 
the ICM, quantitative studies have received comparatively little at- 
tention in hydrodynamical cluster simulations thus far. One reason 
for this is that 3D turbulence is difficult to resolve in any numeri- 
cal scheme, because these always introduce some finite numerical 
viscosity, effectively putting a limit on the Reynolds numbers that 
can still be adequately represented. In the Lagrangian SPH method, 
which has been widely employed for studies of cluster formation, 
an artificial viscosity is used to capture shock s. The original param- 
eterisation of this viscosity ( Monaghan & Gingold.1983) makes the 
scheme comparatively viscous; it smoothes out small-scale velocity 
fluctuations and viscously damps random gas motions well above 
the nominal resolution limit. This hampers the ability of the origi- 
nal SPH to develop fluid turbulence down to the smallest resolved 
scales. 

However, the numerical viscosity of SPH can in principle be 
reduced by using a more sophisticated parameterisation of the ar- 
tificial viscosity. Ideally, the viscosity should only be present in a 
hydrodynamical shock, but otherwise it should be negli gibly small. 
To come closer to this goal. lMorris & M onaghan' fl997) proposed 
a numerical scheme where the artificial viscosity is treated as an in- 
dependent dynamical variable for each particle, with a source term 
triggered by shocks, and an evolutionary term that lets the viscos- 
ity decay in regions away from shocks. In this way, one can hope 
that shocks can still be captured properly, while in the bulk of the 
volume of a simulation, the effective viscosity is lower than in orig- 
inal SPH. We adopt this scheme and implement it in a cosmological 
simulation code. We then apply it to high-resolution simulations of 
galaxy cluster formation, allowing us to examine a better represen- 
tation of turbulent gas motions in SPH simulations of clusters. This 
also shines new light on differences in the results of cosmological 
simulations between different numerical techniques. 

In Section 121 we discuss different ways of implementing the 
artificial viscosity in SPH. We demonstrate in Section|3]the robust- 
ness of our new low-viscosity scheme by applying it to several test 



problems. In Sections El |51 and we introduce our set of cluster 
simulations, the algorithm to detect and measure turbulence, and 
the implications of the presence of turbulence for the structure and 
properties of galaxy clusters. In SectionQ we consider the effects 
of turbulence on the line- width of narrow X-ray metal lines. Fi- 
nally, in Section [S] we apply the results from our new simulations 
to models for the production of radio emission due to turbulent ac- 
celeration processes. We give our conclusions in Section|9l 



2 SIMULATION METHOD 

The smoothed particle hydrodynamics method treats shocks with 
an artificial viscosity, which leads to a broadening of shocks and 
a relatively rapid vo rticity decay. To overcome these problems, 
iMorris & Monaghai] 41997 ) proposed a new parameterisation of 
the artificial viscosity capable of reducing the viscosity in regions 
away from shocks, where it is not needed, while still being able to 
capture strong shocks reliably. We have i mplemented th is method 
in the cosmological SPH code GADGET-2 lSDringei2005ll and de- 
scribe the relevant details in the following. 

In GADGET-2, the viscous force is implemented as 



^ -^^rrijliijS/iWij, 



and the rate of entropy change due to viscosity is 



dt 
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where Ai — {'^ — l)ui / pj~^ h the entropic function of a par- 
ticle of density pi and thermal energy m per unit mass, and 
Wij denotes the arithmetic mean of the two kernels Wij{hi) 
and Wjjjhj). The usual para r neterisation of the artificial viscosity 
iMonaghan & Gingol j FrnsI: lBalsaral[T995.) for an interaction of 
two particles i and j includes terms to mimic a shear and bulk vis- 
cosity. For standard cosmological SPH simulations, it can be writ- 
ten as 



Pij 



(3) 



for rij ■ Vij ^ and Uij = otherwise, i.e. the pair-wise viscosity 
is only non-zero if the particle are approaching each other. Here 



Pij — 



hijVij 



(4) 



Cij is the arithmetic mean of the two sound speeds, pij is the av- 
erage of the densities, hij is the arithmetic mean of the smoothing 
lengths, and rij = Vi — Vj and Vij = Vi — Vj are the inter-particle 
distance and relative velocity, respectively. We have also included 
a viscosity-limiter fij , which is often used to suppress the viscosity 
locally in regions of strong shear flows, as measured by 



(5) 



which can help to avoid spurious angular momentum and vortic- 
ity transport in gas disks (' Steinmet3 ll996). Note however that the 
parameters describing the viscosity (with common choices a — 
0.75 - 1.0, /3 = 2a, ry = O.Ol/i^^, and ai = 0.0001c^//^^ ) stay 
here fixed in time. This then defines the 'original' viscosity scheme 
usually employed in cosmological SPH simulations. We refer to 
runs performed with this viscosity scheme as ovisc simulations. 
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As a variant of the original parameterisation of the artificial 
viscos ity, GADGET-2 can use a formulation proposed by Monasha^ 
based on an analogy with Riemann solutions of compress- 
ible gas dynamics. In this case, /j^ab is defined as 

Vij • Vij 

= (6) 
and one introduces a signal velocity v^^f, for example in the form 



t= 5.0000000 



v-^^ Ci + Cj -Sfiij. 

The resulting viscosity term then changes into 



Pij 



(7) 



(8) 



We have also performed simulations using this signal velocity 
based artificial viscosity and found that it performs well in all test 
problems we examined so far, while in some cases it performed 
slightly better, in particular avoiding post shock oscillations in a 
more robust way. We refer to simulations performed using this 'sig- 
nal velocity' based viscosit y scheme as svisc singulation s. 

The idea proposed bv iMorris & MonaghanI ^1997*) is to give 
every particle its own viscosity parameter a^, which is allowed to 
evolve with time according to 

dai _ ai — t 
"dT ~ r 



(9) 



This causes ai to decay to a minimum value amin with an e-folding 
time T, while the source term Si is meant to make m rapidly 
grow when a particle appro aches a shock. For the decay timescale, 
iMorris & Monaghail il997h proposed to use 



hi I [ci /), 



(10) 



where hi is the smoothing length, a the sound speed and / a free 
(dimensionless) parameter which determines on how many infor- 
mation crossing times the viscosity decays. For an ideal gas and 
a strong shock, this time scale can be related to a length scale 
8 — 0A47/1 (in units of the smoothing length hi) on which the 
viscosity parameter decays behind the shock front. For the source 
term Si, we follow Morris & Monaghan ( 1997) and adopt 



5, = 5*/, max(0,-|(V-^)J), 



(11) 



where (V • v)- denotes the usual SPH estimate of the divergence 
around the particle i. Note that it would in principle be possible 
to use more sophisticated shock detection schemes here, but the 
simple criterion based on the convergence of the flow is already 
working well in most cases. We refer to simulations carried our 
with this 'low' viscosity scheme as Ivisc runs. 

Usually we set 5** =0.75 and choose / = 1. We also restrict 
ai to be in the range amin = 0.01 and amax = 0.75. Choos- 
ing amin > has the advantage, that possible noise which might 
be present in the velocity representation by the particles on scales 
below the smoothing length still will get damped with time. In- 
creasing S* can give a faster response of the artificial viscosity to 
the shock switch without inducing higher viscosity than necessary 
elsewhere. We also note that we replace a in equation |3] (and equa- 
tion |8] respectively) by the arithmetic mean aij of two interacting 
particles. Depending on the problem, we initialise ai at the start 
of a simulation either with ttmin or ttmax, depending on whether 
or not there are already shocks present in the initial conditions, re- 
spectively. 

While the approach to reduce numerical viscosity with a time- 
variable ai works well with both basic parameterisations of the ar- 
tificial viscosity, most of our cosmological simulations were carried 




Figure 1. A standard shock tube problem computed with the 

low-viscosity scheme with an individual, time-dependent viscosity. From 
top to bottom, we show current value of the strength of the artificial vis- 
cosity ai, density, velocity, pressure, and internal energy, averaged for bins 
with spacing equal to the SPH smoothing length for particles in the low 
density region. The analytic solution of the problem for the time t = 5.0 is 
shown as a solid line. 

out with the 'original' parameterisation because the signal velocity 
variant became available in GADGET-2 only recently. 



3 TEST SIMULATIONS 

To verify that the low-viscosity formulation of SPH with its time- 
dependent artificial viscosity is still able to correctly capture strong 
shocks, we computed a number of test problems. We here report on 
a standard shock tube test, and a spherical collapse test, which both 
have direct relevance for the cosmological formation of halos. As 
a more advanced test for the ability of the code to follow vorticity 
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velocity 



entropy 





Figure 2. Profiles of velocity (left column), pressure (middle left column) entropy (middle right column) and viscosity constant a (right column) for the 
spherical collapse test at 2 different times (from top to bottom). The thin line marks the analytic solution, diamonds give the result obtained by the new SPH 
formulation for the time dependent viscosity. The thick line is the analytic solution adaptively smoothed with the SPH kernel, using the smoothing length of 
the particles in each bin. The lengths of the horizontal lines plotted at each data point correspond to the smoothing lengths of the SPH particles at this position. 



generation, we investigate the problem of a strong shock striking an 
overdense cloud in a background medium. This test can also give 
hints whether turbulence is less strongly suppressed in the low- 
viscosity treatment of SPH than in the original formulation. 



3.1 Shock-Tube Test 

First, we computed a shock tube proble m, which is a common test 
for hydrodynamical numerical schemes llSodI For definite- 

ness, a tube is divided into two halves, having density pi = 1 and 
pressure pi = 1 on the left side, and p2 = 0.125 andp2 = 0.1 on 
the right side, respectively. Like in Sod (1978), we assume an ideal 
gas with 7 = 1.4. To avoid oscillations in the post shock region 
(note that a shock is present in the initial conditions) we initialise 
the viscosity values of the particles with amax = 0.75. We com- 
pute the test in 3D and make use of periodic boundary conditions. 
The initial particle grid is 5x5x50 on one half, and 10x10x100 on 
the other half, implying an equal particle mass for both sides. 

In Figure [I] we show the state of the system at simulation 
time t = 5 in terms of density, velocity, and internal energy with 
a binning which corresponds to the smoothing length for particles 
in the low density region. We also include the analytic expecta- 
tion for comparison. In addition, we plot the values of the artificial 
viscosity parameter of the particles. Clearly visible is that the vis- 
cosity is close to ttmin everywhere, except in the region close to the 
shock. One can also see how the viscosity builds up to its maximum 
value in the pre- shock region and decays in the post shock region. 
We note that the final post- shock state of the gas agrees well with 
the theoretical expectation, and is indistinguishable from the case 
where the original viscosity parameterisation is used. 



3.2 Self-similar spherical collapse 

A test arguably more relevant for cosmological structure formation 
is the self- similar, spherical collapse of a point perturbati on in a ho- 
mogeneous expanding background ('Bert schingeilll985l) . This test 
is difficult for grid and SPH codes alike. The gas cloud collapses 
self similarly, forming a very strong shock (with formally has in- 
finite Mach number) at the accretion surface. While grid codes 
with shock capturing schemes can usually recover the sharp shock 
surface very well, the large dynamic range in the post shock re- 
gion with its singular density cusp, as well as the strict spherical 
symmetry of the problem, are challenging for mesh codes. On the 
other hand, Lagrangian SPH codes tend to have fewer problems 
with the central structure of the post- shock cloud, but they broaden 
the shock surface substantially, and typically show appreciable pre- 
shock entropy injection as result of the artificial viscosity. 

We have computed the self- similar collapse test and compared 
the results for the new viscosity parameterisation with the ana- 
lytic expectation. The very strong spherical shock of this problem 
is a particularly interesting test, because we can here test whether 
the low-viscosity formulation is still able to capture the strongest 
shocks possible. 

In Figure |2| we show the structure of the shock at 2 consecu- 
tive times, scaled to the self-similar variables. In general, the SPH 
result recovers the analytic solution for the post- shock state very 
well, especially when the entropy profile is considered. However, 
the shock is substantially broadened, and some pre-heating in front 
of the shock is clearly visible. In the velocity field, some weak 
post-shock oscillations are noticable. We have also indicated the 
smoothing lengths of the SPH particles as horizontal error bars for 
each of the data points (the points at which the SPH kernel falls to 
zero is reached at twice this length). For comparison, we addition- 
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Figure 3. Time evolution of the interaction of a strong shock wave with an overdense cloud. We show the projected gas density and compare simulations 
carried out with original SPH (left) and the low-viscosity formulation (right). The incident shock wave has Mach number 10, and the cloud is initially at 
pressure equilibrium in the ambient medium and has overdensity 5. 



ally over-plotted the analytic solution adaptively smoothed with the 
SPH kernel size at each bin. 

The panels of the right column in Figure |2| show the profile 
of the viscosity parameter, which was set to amin at the beginning 
of the simulation, as the initial conditions do not contain a shock. 
The viscosity parameter builds up immediately after starting the 
simulation as the strong shock forms. Later one can see how the 
viscosity parameter begins to decay towards amin in the inner part, 
how it builds up to amax towards the shock surface, and how a 
characteristic profile develops as the shock moves outward. In the 
post- shock region an intermediate viscosity values is maintained 



for some time due to some non-radial motions of gas particles in 
this region. 



3.3 Schock-cloud interaction 

To verify that the low-viscosity scheme also works in more com- 
plex hydrodynamical situations, we simulate a test problem where a 
strong shock strikes a gas cloud embedded at pressure equilibrium 
in a lower density environment. A recent discussion of this setup 
can be found in Poludnenko et al. (2002) and references therein. 
SPH is able to reproduce the main features expected in this test 
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non radiative 

Table 1. Main characteristics of the non radiative galaxy cluster simulations. Column 1: identification label. Columns 2 and 3: mass of the dark matter (Mdm) 
and gas (Mgas) components inside the virial radius. Column 4: virial radius R^. Column 5: X-ray luminosity inside the virial radius Lx- Columns 6 and 7: 
mass-weighted temperature (Tmw) and spectroscopic like temperature (Tsl)- 



Simulations 


MDM(/i" 








R^{h- 


1 kpc) 




ergs ■*■) 


TMw(keV) 




(keV) 




svisc 


Ivisc 


svisc 


Ivisc 


svisc 


Ivisc 


svisc 


Ivisc 


svisc 


Ivisc 


svisc 


Ivisc 


gl 


14.5 


14.5 


17.5 


17.0 


2360 


2355 


47.1 


21.3 


7.2 


7.1 


5.8 


5.6 


g8 


22.6 


22.4 


19.8 


19.8 


2712 


2705 


63.1 


32.1 


9.3 


9.1 


6.2 


5.7 


g51 


13.0 


13.0 


11.5 


11.5 


2255 


2251 


30.8 


17.9 


6.4 


6.3 


4.6 


4.7 


g72 


13.5 


13.4 


12.0 


11.9 


2286 


2280 


18.3 


14.1 


5.8 


5.8 


4.0 


4.0 


g676 


1.1 


1.0 


0.95 


0.91 


983 


972 


3.2 


1.4 


1.3 


1.3 


1.6 


1.5 


g914 


1.2 


1.0 


1.07 


0.91 


1023 


971 


4.2 


1.7 


1.4 


1.3 


1.6 


1.7 


gl542 


1.1 


1.0 


0.95 


0.90 


982 


967 


3.0 


1.4 


1.3 


1.2 


1.4 


1.5 


g3344 


1.1 


1.1 


1.00 


0.96 


1002 


993 


2.2 


1.4 


1.4 


1.3 


1.4 


1.5 


g6212 


1.1 


1.1 


1.00 


1.01 


1000 


1006 


3.0 


1.5 


1.3 


1.3 


1.6 


1.7 



problem reasonably well, like reverse and reflected shocks, back- 
flow, primary and secondary Mach stems, primary and secondary 
vortices, etc. (see Sori ngell 120051) . Our purpose here is to check 
whether the new scheme for a time-variable viscosity performs at 
least equally well as the original approach. 

In Figure |3l we compare the time evolution of the projected 
gas density for the original viscosity scheme (left hand side) with 
the new low- viscosity scheme (right hand side). Overall, we find 
that the new scheme produces quite similar results as the origi- 
nal method. But there are also a number of details where the low- 
viscosity scheme appears to work better. One is the external reverse 
bow shock which is resolved more sharply with the new scheme 
compared to the original one. This is consistent with our findings 
from the previous tests, where we could also notice that shocks 
tend to be resolved somewhat sharper using the new scheme. We 
also note that instabilities along shear flows (e.g. the forming vor- 
texes or the back-flow) are appearing at an earlier time, as expected 
if the viscosity of the numerical scheme is lower. This should help 
to resolve turbulence better. 

In summary, the low-viscosity scheme appears to work with- 
out problems even in complex situations involving multiple shocks 
and vorticity generation, while it is still able to keep the advan- 
tage of a reduced viscosity in regions away from shocks. We can 
therefore expect this scheme to also work well in a proper environ- 
ment of cosmological structure formation, and simulations should 
be able to benefit from the reduced viscosity characteristics of the 
scheme. 



sured that all of our clusters are free from contaminating boundary 
effects out to at least 3 - 5 virial radii. Gas was introduced in the 
high-resolution region by splitting each parent particle into a gas 
and a DM particle. The final mass-resolution of these simulations 
was rriBM = 1.13 x 10^ H'^Mq and mgas = 1.7 x 10^ H'^Mq 
for dark matter and gas within the high-resolution region, respec- 
tively. The clusters were hence resolved with between 2 x 10^ 
and 4 x 10^ particles, depending on their final mass. For details 
on their properties see Table [I] The gravitational softening length 
was e = 5.0 /i~^kpc (Plummer-equivalent), kept fixed in physical 
units at low redshift and switched to constant comoving softening 
of e = 30.0/i~^kpc at z ^ 5. Additionally we re-simulated one 
of the smaller cluster (g676) with 6 times more particles (HR), de- 
creasing the softening by a factor of two to e = 2.5 h~^kpc. 

We computed three sets of simulations using non radiative gas 
dynamics, where each cluster was simulated three times with dif- 
ferent prescriptions for the artificial viscosity. In our first set, we 
used the original formulation of artificial viscosity within SPH. In 
the second set, we used the parametrisation based on signal veloc- 
ity, but with a fixed coefficient for the viscosity. Finally, in our third 
set, we employed the time dependent viscosity scheme, which we 
expect to lead to lower residual numerical viscosity. Our simula- 
tions were all c arried out with an extended version of GADGET-2 
(*SDri ngei2005h . a new version of the parallel TreeSPH simulation 
code GADGET fSDringel et al. 2001). We note that the formulation 
of SPH used in this code follows t he 'ent ropy-conserving' method 
proposed bv Soringel & Hernaui's] i2002h . 



4 COSMOLOGICAL CLUSTER SIMULATIONS 

We have performed high-resolution hydrodynamical simulations of 
the formation of 9 galaxy clusters. The clusters span a mass-range 
from 10^^ h~^MQ to 2.3 x lO^^h'^Mc^ and have originally been 
selected from a DM-only simulation ( Yoshida et al. 2001) with 
box-size 479/i"^Mpc of a flat ACDM model with Qo = 0.3, 
h = 0.7, as = 0.9 and = 0.04. Using the 'Zoomed Ini- 
tial Conditions' (ZIC) technique jTormen et alJll997n . we then re- 
simulated the clusters with higher mass and force resolution by 
populating their Lagrangian regions in the initial conditions with 
more particles, adding additional small-scale power appropriately. 
The selection of the initial region was carried out with an iterative 
process, involving several low resolution DM-only resimulations to 
optimise the simulated volume. The iterative cleaning process en- 



5 IDENTIFYING TURBULENCE 

In the idealized case of homegeneus and isotropic turbulence, the 
autocorrelation function of the velocity field of the fluid should not 
depend on the position (homogeneity) and it should only depend 
on the magnitude of the distance r between points (isotropy). The 
tens or of the correlati on function of the velocities is thus given by 
(e.g.lChoudhurilElli): 

R^j{r) = {v^{x)vj{x + r)) (12) 

where x is the position of a fluid particle. The 3D power spectral 
density of the turbulent field is given by (e.g. C houdhuri 1998) : 

^ii(fe) = tAi [ fii, (r)exp(ifer)dr. (13) 
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Figure 4. Mean local velocity dispersion for the central 500^ kpc^ box as a 
function of the resolution adopted for the TSC-smoothing of the local mean 
field. Results are plotted for a low viscosity simulation. 



The energy spectrum, E{k), associated with the fluctuations of the 
velocity field is related to the diagonal parts of both the tensor of 
the correlation function, and that of the power spectral density. This 
energy spectrum is given by (e.g. .ChoudhurL 1 99 : 



E{k) =27Tk^^u{k), 

and the total turbulent energy per unit mass is 



'^^turb 



E{k) dk, 



(14) 



(15) 



where the summation convention over equal indices is adopted. 

The real case of the intracluster medium is however much 
more complex, in particular not homogeneous and isotropic. The 
gravitational field induces density and temperature gradients in the 
ICM, and the continuous infall of substructures drives bulk mo- 
tions through the ICM. These effects break both homogeneity and 
isotropy at some level, at least on the scale of the cluster, and thus 
demand a more complicated formalism to appropriately charac- 
terise the turbulent field. It is not the aim of the present paper to 
solve this problem completely. Instead we focus on a zero-order 
description of the energy stored in turbulence in the simulated 
boxes, and for this purpose the basic formalism described below 
should be sufficient. 

A crucial issue in describing turbulent fields in the ICM is the 
distinction between large-scale coherent velocity field and small- 
scale 'random' motions. Unfortunately, the definition of a suit- 
able mean velocity field is not unambiguous because the involved 
scale of averaging introduces a certain degree of arbitrariness. Per- 
haps the simplest possible procedure is to take the mean velocity 
computed for the cluster volume (calculated, for example, within 
a sphere of radius Rvir) as the coherent velocity field, and then 
to define the turbulent velocity component as a residual to this 
velocity. This simple approach (hereafter stand ard approach) has 
been widely employed in previous works fe.g. JNorman & BrvarJ 



Il999l: ISunvaev et alJl2003h . and led to the identification of ICM 
turbulence in these studies. However, an obvious problem with this 
method is that this global subtraction fail to distinguish a pure lam- 
inar bulk flow from a turbulent pattern of motion. Note that such a 
large scale laminar flows are quite common in cosmological sim- 
ulations, where the growth of clusters causes frequent infalls and 
accretions of sub-halos. This infall of substructures is presumably 
one of the primary injection mechanisms of ICM turbulence. 

To avoid this problem, a mean velocity field smoothed on 
scales smaller than the whole box can be used, and then the field of 
velocity fluctuations is defined by subtracting this mean-local ve- 
locity, v{x), from the individual velocities Vi of each gas particle. 
We note that if the smoothing scale is chosen too small, one may 
risk loosing large eddies in the system if they are present, but at 
least this procedure does not overestimate the level of turbulence. 

Following this second approach (hereafter local-velocity ap- 
proach), we construct a mean local velocity field v{x) on a uni- 
form mesh by assigning the individual particles to a mesh with a 
Triangular Shape Cloud (TSC) window function. The mesh covers 
a region of 1.0 comoving Mpc on a side and typically has between 
8^ and 64^ cells, which is coarse enough to avoid undersampling 
effects. The equivalent width of the TSC kernel is approximatively 
3 grid cells in each dimension, corresponding to a smoothing scale 
of ^ 360 — 45 kpc, respectively. As our analysis is restricted only 
to the highest density region in the clusters, the scale for the TSC- 
smoothing is always larger than the SPH smoothing lengths for the 
gas particles, which typically span the range 7.5 — 15 /i~^kpc in 
the box we consider. 

We then evaluate the local velocity dispersion at the position 
X of each mesh cell over all particles a in the cell by: 



j{x) {[Va,i - Vi{Q 



(^)l)c 



(16) 



where i and j are the indices for the three spatial coordinates, and 
Oceii denotes the average over particles within each cell. 
The diagonal part of the tensor of the correlation function of the 
field of velocity fluctuations at r = in the simulated box can then 
be approximated by 

7?«(r = 0)~(a?,(a:))g^^. (17) 

Based on Equation I15i . we can then estimate the energy density of 
the turbulence in real space as 



p{x) 



E{k) dkr.-p{x) 



(18) 



in the local-velocity and standard case, respectively. Here p{x) is 
the gas density within the cells. 

The subtraction of a local velocity from the velocity distribu- 
tion of the particles is expected to efficiently filter out the contribu- 
tion from laminar bulk-flows with a scale ^ 3 times the size of the 
cells used in the TSC smoothing. However, a large-scale turbulent 
velocity field component, if it exists, would also be suppressed, so 
that this procedure can be expected to reduce the turbulent velocity 
field to a certain degree. As shown in Figure|4] this depends on the 
resolution of the mesh used in the TSC assignment. Fig. |4l shows 
that the increase of the turbulent velocity dispersion with the cell 
size is not dramatic for cell sizes larger than 100 kpc. We find that 
(Vazza et al., in prep.) a TSC smoothing with larger cell sizes would 
not efficiently filter out contributions from laminar bulk-motions. It 
can be tentatively concluded that the local velocity approach with 
a smoothing with 16^ - 32^ cells in the central (1.0 Mpc)^ vol- 
ume catches the bulk of the turbulent velocity field in the simulated 
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Figure 5. Gas velocity field in a slice through the central Mpc of a cluster simulation g72 after subtracting the global mean bulk velocity of the cluster. The 
panels on the left is for a run with the original viscosity of SPH while the panels on the right shows the result for the low viscosity scheme. The underlying 
colour maps represent the turbulent kinetic energy content of particles, inferred using the local velocity method (upper row) or the standard velocity method 
(lower row). For the local velocity method a conservative 64^ grid is used in the TSC smoothing. The cluster centre is just below the lower-left corner of the 
images. The vertical lines in the upper row show where the 1-dimensional profile for the simulated radio-emission of Fig.[l9|are taken. 



box. Therefore, if not specified otherwise, all the numerical quan- 
tities given in the following are obtained using a TSC-assignment 
procedure based on 32^ cells. A more detailed discussion of this 
method and tests of the parameters involved is reported elsewhere 
(Vazza et al., in prep.). 

Figures |5] and |6l give examples of the turbulent velocity field 
calculated with both the standard and local velocity methods, show- 
ing the same galaxy cluster in both cases, but in one case simulated 
with the signal-velocity variant of the viscosity, and in the other 



with the new time-dependent low-viscosity scheme. Note that we 
here selected a situation where a large (ca. 500 kpc long) lami- 
nar flow pattern can be easily identified close to the centre of one 
of our simulated clusters {g72). When the mean cluster velocity 
field is subtracted as in Figure 13 large residual bulk flow patterns 
remain visible, caused by a substructure moving through the clus- 
ter atmosphere. We colour-coded the turbulent kinetic energy of 
particles, Et{x) ~ 1/2 p(x)av(xy, after subtracting the local 
mean velocity field (here smoothed onto a 64^ mesh) for the up- 
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Figure 6. Same slice of the Gas velocity field as in figure|5|of cluster g72 after subtracting the local mean velocity of the cluster. The panel on the left is for a 
run with the original viscosity of SPH while the panel on the right shows the result for the low viscosity scheme. 



per panels and aster subtracting the global mean bulk velocity of 
the cluster for the lower panels. One can see that fluid instabili- 
ties of Kelvin-Helmholtz type are growing along the interfaces of 
the large laminar flow pattern, visible in the upper left panel. As 
expected, the strength of this turbulent velocity field is consider- 
ably larger in the simulation obtained with the new low-viscosity 
scheme, providing evidence that such instabilities are less strongly 
damped in this scheme. This can also be seen by the longer flow 
field lines in Figure |6l Figures |5] also visually confirms the differ- 
ences in the two approaches of filtering the velocity field. Whereas 
the local-velocity approach highlights the energy within the veloc- 
ity structure along boundary layers, the energy within the large, 
bulk motions are preferentially selected when only subtraction the 
global mean bulk velocity. 

The total cumulative kinetic energy in the random gas motions 
inside our mesh (cantered on the cluster centre) reaches 5%-30% of 
the thermal energy for the simulations using the new, low-viscosity 
scheme, whereas it stays at much lower levels (^2%-l0%) when 
the signal velocity parameterisation of the viscosity is used. If the 
original viscosity scheme is used, it is typically at even lower values 

In general, we find that more massive clusters tend to have a 
higher fraction of turbulent energy content. However, given that our 
simulations have fixed mass resolution, this trend could in princi- 
ple simply reflect a numerical resolution effect. In order to get fur- 
ther information on this, we have re-simulated one of the smaller 
clusters (g676) with 6 times better mass resolution using the signal 
velocity parameterisation of the viscosity. Ai z — 0, this cluster is 
then resolved by nearly as many particles as the massive clusters 
simulated with our standard resolution. We find that for this high- 
resolution simulation the level of turbulence 3%) is increased 
compared with the normal resolution 2%), but it stays less to 
what we found for the low-viscosity scheme at our normal resolu- 
tion 5%). This confirms two expectations. First, the low viscos- 
ity scheme effectively increases the resolution on which SPH simu- 



lations can resolve small-scale velocity structure, which otherwise 
gets already suppressed on scales larger than the smoothing length 
by spurious viscous damping effects due to the artificial viscosity. 
Second, the amount of turbulence in the high resolution version of 
g676 is still less than what we find with the same viscosity imple- 
mentation in the larger systems, and even much smaller than what 
we find with the low- viscosity scheme in the large clusters. This 
tentatively suggests that the trend of a mass-dependence of the im- 
portance of turbulence is not caused by numerical effects. Note that 
with a fixed physical scale of 1 Mpc we are sampling different frac- 
tions of i?vir in clusters of different masses. However, if, in case of 
the less massive system, we restrict the sampling relative to i?vir to 
measure within comparable volumes, the fraction of turbulent en- 
ergy content found in the small cluster increase roughly by a factor 
of two. Thereby we still find a significant trend with mass when 
measuring turbulence within a fixed fraction of Rvir- Although it 
should be mentioned, that unless the dissipation of turbulence on 
small scales will me modeled correctly in a physical granted way, 
the different formation time scales of systems with different masses 
can potentially also contribute to such a trend. 

In order to verify that our method for measuring the local ve- 
locity dispersion gives reasonable values. Figure shows a radial 
profile of the volume-weighted, relative difference between ther- 
mal pressure for the signal velocity based and low-viscosity run. 
Here we used the an average over the three massive clusters (gl,g51 
and g72) which have comparable masses. The solid line shows the 
relative difference in radial shells and indicates that the turbulent 
pressure support can reach even up to 50% in the central part and 
drops to at approximately 0.2 Rvir. The dashed line shows the 
cumulative difference, which over the total cluster volume con- 
tributes between 2% and 5% to the total pressure. The diamonds 
mark the measurement inferred from the local velocity dispersion 
within centred boxes of various sizes. We also calculate the dif- 
ference between the signal velocity based and low-viscosity runs 
using the mean values over the three clusters. Qualitatively, there 
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Figure 7. Radial profile of the relative thermal pressure difference averaged 
over three nearly equally massive clusters (gl,g51 and g72), comparing the 
signal velocity based and low-viscosity runs (lines). The dashed line is the 
cumulative difference, whereas the solid line marks the profile in radial 
shells. The diamonds mark the difference in the turbulent energy support we 
inferred from the local velocity dispersion within several concentric cubes 
of different sizes (/cube = 2r) for the same runs. This should be compared 
with the dashed line. The inlay shows the absolute value inferred from the 
local velocity dispersion from the different viscosity parameterisations, re- 
spectively. 



is good agreement of results obtained with this approach with the 
cumulative curve. This confirms that our method to infer the turbu- 
lent energy content from the local velocity dispersion of the gas is 
meaningful. Note that the temperature which is used to calculate the 
pressure is determined by strong shock heating. As different resim- 
ulations of the same object can lead to small (but in this context 
non-negligible) timing differences, this can introduce sizable vari- 
ations in the calculated pressure, especially during merging events. 
We verified that these differences for individual clusters are sig- 
nificantly larger than the differences between the cumulative curve 
(dashed line) and the data points from the local velocity dispersion 
(diamonds). Therefore we can only say that the two methods agree 
well within their uncertainties. 

Finally, the inlay of Figure gives the absolute contribution 
from the low-viscosity, the original viscosity in its two variants us- 
ing the local velocity dispersion respectively. It seems that using the 
signal based viscosity in general leads already to more turbulence 
than the "old" original viscosity, but the time-dependent treatment 
of the viscosity works even more efficiently. 

Although we are using a formalism which is suitable only 
for isotropic and homogeneous turbulence, the study of the turbu- 
lent energy spectrum may provide some useful insight. In the lo- 
cal mean velocity approach, we can obtain the diagonal part of the 
turbulent energy spectrum using Equation IT3l . with Ru approxi- 
mated as 



Rii(r) = {[Va,i - Vi{Xa)] [Vb,i " Vi{Xb)])^^ 



(19) 



Figure 8. The energy spectra of the standard velocity fluctuations (upper 
curves) and of the local velocity fluctuations (lower curves) of gas parti- 
cles in the central 500^ kpc^ region of a cluster simulated with the original 
recipe for the artificial viscosity, with signal-velocity and with the low vis- 
cosity implementation. Additionally a Kolmogorov slope (dot-dot-dashed) 
is drawn for comparison. 



where v{xa) is the TSC-mean velocity of the cell which contains 
the point Xa, and the average is over all pairs (a, b) in the box 
with a certain distance r. In the standard approach, we would here 
subtract the centre-of-mass velocity of the cluster instead. 

A major problem for estimating the correlation functions 
Rii(r) in this way, and with the energy spectrum calculated from 
SPH simulations (and in general from adaptive resolution ap- 
proaches), is given by the non-uniform sampling of the point 
masses in the simulated box. To reduce this problem we focus on 
regions corresponding to the cores of galaxy clusters. Here the re- 
quirement of isotropic and homogeneous turbulence is hopefully 
better full filled. Also the density profile is relatively flat such that 
the sampling with gas particles is not too non-uniform. In addition, 
we estimate the correlation function as an average of dozens of 
Monte-Carlo extractions of gas particles from the simulated out- 
put, where we picked one particle from each of the (15.6 kpc)^ 
cells in order to have a uniform, volume- weighted set of particles. 

In Figure |8] we show examples for the energy spectra we ob- 
tained for the two approaches by Fourier-transforming the mea- 
surements fov Rii(r). The energy spectra for the two methods for 
treating the mean velocities are reasonably similar in shape, but 
the spectrum calculated with the local mean velocity has a lower 
normalisation, independent of resolution. This is expected because 
the TSC smoothing filters out contributions from laminar motions, 
and may also damp the turbulent field at some level. Both spectra 
show a slope nearly as steep as a Kolmogorov spectrum (which has 
E(k) oc k~^^^) at intermediate scales, but exhibit a significant flat- 
tening at smaller scales (i.e. large k). The flattening at small scales 
could be caused by numerical effects inherent in the SPH tech- 
nique, where an efficient energy transfer to small-scale turbulent 
cells on scales approaching the numerical resolution is prevented. 
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and thus a complete cascade cannot develop. Additional, the lack 
of numerical viscosity in the low-viscosity scheme can in princi- 
ple lead to an increase of the noise level within the velocity field 
representation by the SPH particles on scales below the smoothing 
length. Such noise in general could contribute to the flattening at 
small scales. It is however not clear how to separate this noise com- 
ponent from a turbulent cascade reaching scales similar or below 
the smoothing length. Therefore one focus of future investigations 
has clearly be towards this issue. It is however important to note 
that the largest turbulent energy content (expecially at small scales) 
is always found in the clusters simulated with the low-viscosity 
scheme. This is particularly apparent in the energy spectra when 
the local velocity approach is used and suggests that the energy 
spectrum obtained with the standard approach is significantly af- 
fected by laminar bulk-flows, which are not sensitive to a change 
in parameterisation of the artifical viscosity. 



6 CLUSTER PROPERTIES 

Different levels of small-scale random gas motions within the ICM 
have only mild effects on global properties of clusters like mass or 
temperature, as evidenced by the measurements in Table [l] How- 
ever, additional kinetic energy in turbulent motions changes the 
central density and entropy structure, which in turn has a sizable 
effect on the X-ray luminosity. We investigate the resulting changes 
in cluster scaling relations and radial profiles in more detail within 
this section. 

6.1 Maps 

The presence of a significant turbulent energy component in the 
intra-cluster medium manifests itself in a modification of the bal- 
ance between the gravitational potential and the gas pressure. There 
are in fact observational reports that claim to have detected such 
fluctuations i n pressure maps derived from X-ray observations 
iSchuecker et al. 2004). We here calculate artificial pressure Part 
maps for our simulation s, based on surf ace brightness maps (Lx) 
and spectroscopic-like iMazzotta et alJ |2004) temperature (Tsi) 
maps. This allows artificial pressure maps to be estimated as 

Part^nTsi, (20) 

where we defined 

n= {Lx/VtTi)^ . (21) 

Figures l9l and fTol show a comparison of a number of cluster 
maps produced using an unsharp-mask technique of the form 

Image^^^harpmask = I^age - Smoothed (Image, a), (22) 

where a Gaussian smoothing with FWHM of a = 200 kpc 
was applied. We analyse maps of the X-ray surface brightness, 
spectroscopic-like temperature, 'true' pressure maps (e.g. based on 
Compton y) and artificial pressure maps constructed as described 
above. All maps show the central 2 Mpc of the cluster run gl, sim- 
ulated with the low-viscosity scheme (right panels) compared with 
the original SPH scheme (left panels). 

Disregarding the large contribution by substructure in all the 
X-ray related maps (therefore also in the artificial pressure map), all 
types of maps show clear signs of turbulence. It is noticeable in both 
runs, but it has a much larger extent and intensity in the the low- 
viscosity run. Note in particular the turbulent motions (appearing 
as lumpiness in the unsharp-mask images) in the wake of infalling 



substructures, and the earlier break-up of fluid interfaces when the 
new, reduced viscosity scheme is used. 

Pressure maps (and therefore SZ maps) are arguably the most 
promising type of map when searching for observational imprints 
of turbulence. Apart from reflecting the large scale shape of the 
cluster they are known to be relatively featureless, because most of 
the substructures in clusters are in pressure equilibrium with their 
local surroundings, making them in principle invisible in pressure 
maps. On the other hand, the contribution of the turbulent motion to 
the local pressure balance can be expected to leave visible fluctua- 
tions in the thermal pressure map. This can indeed be seen nicely in 
the pressure (e.g. SZ) maps in Figure[lO| Note that the amplitudes 
of the turbulent fluctuations in the case of the low-viscosity run 
are larger and also spatially more extended in the core of the clus- 
ter. Artificial pressure maps constructed from the X-ray observables 
still show such fluctuations, but they are partially swamped by the 
signatures of the infalling substructures, making it difficult to quan- 
tify the amount of turbulence present in clusters using such X-ray 
based artificial pressure maps. 

The small displacements seen in the substructure positions 
between the two runs are due to small differences in their orbital 
phases. Besides the general problem to precisely synchronise clus- 
ter simulations w ith different dynamical processes involved , it is 
well known fe.g. iTormen et al J 120041: IPuchwein et alJl2005h that 
the interaction of the gas with its environment can significantly 
change the orbits of infalling substructure. The different efficien- 
cies in stripping the gas from the infalling substructure in the sim- 
ulations with different viscosity prescription can therefore lead to 
small differences in the timing and orbits between the two simula- 
tions. 



6.2 Scaling Relations 

In Figure [n] we compare the mass-weighted temperature of our 
galaxy clusters for simulations with the original viscosity and 
for runs with the low-viscosity scheme. There are no significant 
changes. Comparing the X-ray luminosity instead, we find that it 
drops significantly by a factor of ^ 2 for clusters simulated with 
the low- viscosity scheme, as shown in Figure [l2l This is quite in- 
teresting in the context of the long-standing problem of trying to 
reproduce the observed cluster scaling relations in simulations. In 
particular, since non radiative cluster simulations tend to produce 
an excess of X-ray luminosity, this effect would help. However, 
one has to keep in mind that the inclusion of additional physical 
processes like radiative cooling and feedback from star formation 
can have an even larger impact on the cluster luminosity, depend- 
ing on cluster mass, so a definite assessment of the scaling relation 
issue has to await simulations that also include these effects. 



6.3 Radial profiles 

The presence of turbulence manifests itself in an increase of the ve- 
locity dispersion of the cluster gas, especially towards the centre, 
while the dark matter velocity dispersion should be unaffected. In 
Figure[n| we compare the velocity dispersion of gas and dark mat- 
ter, averaged over the low- and high-mass clusters in our set. As ex- 
pected, the velocity dispersion of the dark matter does not change 
in the low-viscosity simulations, where a larger degree of turbu- 
lence is present in the ICM. On the other hand, the central velocity 
dispersion of the gas increases, reaching of order of 400 km for 
our massive clusters. As the gas is in pressure equilibrium with the 
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Figure 9. Unsharp mask images of X-ray maps for one of the massive clusters (g7), comparing runs with the low-viscosity scheme (right panels) with the 
original SPH scheme (left panels). The top row gives maps of surface brightness, while the bottom one compares maps of the 'spectroscopic like' temperature, 
both within 2 Mpc centred on the cluster. We can see evidence for an increased level of turbulent motions behind the infalling substructures, and a break-up of 
fluid interfaces for the reduced viscosity scheme is clearly visible. Also, there is a general increase of turbulence (appearing as lumpiness) towards the centre. 
However, the most prominent signals in the map stem from the higher density or different temperature of substructures relative to their surrounding, or from 
shocks and contact discontinuities. For this reason, turbulence can be better identified in pressure maps (see FigurelTol. 



unchanged gravitational potential, the hydrodynamic gas pressure 
will be correspondingly lower in the centre due to the presence of 
these random gas motions. 

In Figure [21 we show the mean cluster temperature profiles, 
which only shows a very mild trend of increasing temperature in the 
central part of clusters when using the new, low-viscosity scheme. 
However, the central gas density drops significantly in the low- 
viscosity scheme, as shown in FigurelTH This change in density is 
restricted to inner parts of the cluster, roughly to within ^XRvir, 



which may be taken as an indication of the size of the region where 
turbulent motions are important. 

Quite interestingly, the presence of turbulence also changes 
the entropy profiles. In Figure [l6l we show the radial entropy pro- 
files of our clusters, which in the case of the low-viscosity scheme 
exhibit an elevated level of entropy in the core, with a flattening 
similar to that inferred from X-ray observations. It is remarkable 
that this central increase of entropy occurs despite the fact that the 
source of entropy generation, the artificial viscosity, is in princi- 
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Figure 10. Unsharp mask images of pressure maps of one of the massive clusters (gl), comparing runs with the low-viscosity scheme (right panels) with 
the original SPH scheme (left panels). We also compare different methods for determining the pressure maps. The panels of the top row show Compton-y 
maps (which can be associated with projected, thermal pressure maps), whereas the maps in the bottom row are pressure maps derived based on X-ray surface 
brightness and spectroscopic temperature maps, see equation |20|and|2] Both kinds of maps show an increase of structure (lumpiness) for the simulation which 
uses the reduced viscosity scheme (right panels) when compared with the original SPH viscosity scheme (left panels). The maps based on X-ray observables 
show a larger degree of lumpiness due to the gas around substructures, especially in the vicinity of infalling subgroups. 



pie less efficient in the low-viscosity scheme. There are two main 
possibilities that could explain this result. Either the low-viscosity 
scheme allows shocks to penetrate deeper into the core of the clus- 
ter and its progenitors such that more efficient entropy production 
in shocks occurs there, or alternatively, the reduced numerical vis- 
cosity changes the mixing processes of infalling material, allowing 
higher entropy material that falls in late to end up near the cluster 
centre. 

In order to investigate a possible change of the accretion be- 



haviour, we traced back to high redshift all particles that end up at 
z = within 5% of i?vir of the cluster centre. We find that most 
of the central material is located in the centres of progenitor halos 
at higher redshift, which is a well known result. However, in the 
simulations with the time dependend, low-viscosity scheme, there 
is a clear increase of the number of particles which are not asso- 
ciated with the core of a halo at higher redshift. We illustrate this 
with the histograms shown in Figure [TTl which gives the distribu- 
tion of the distance to the nearest halo in units of i^vir of the halo. 
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Figure 11. Comparison of the virial temperature of the 9 clusters when 
different parameterisations of the viscosity are employed. The solid line 
marks the one-to-one correspondence. 



All particles at distances larger than 1 are not associated with any 
halo at corresponding epoch. Compared to the low entropy material 
that is already bound in a dense core at this epoch, this diffuse gas 
is brought to much higher entropy by shocks. When it is later ac- 
creted onto the cluster and mixed into the core, it can then raise the 
entropy level observed there. We note that Eulerian hydrodynamics 
simulations also show a flattening of the entropy profile. While the 
exact degree to which numerical and physical (turbulent) mixing 
contribute to producing this result is still a matter of debate, it is 
intriguing that a larger level of turbulence in the SPH calculations 
substantially alleviates the discrepancies in the results otherwise 
obtained with the two techniques iFrenk et aLJ999. : ,^0'Shea et alJ 
|20 



Figure 12. Comparison of the bolometric luminosity of the 9 clusters when 
different parameterisations of the viscosity are employed. The solid line 
marks the one-to-one correspondence. It is evident that clusters with a larger 
degree of turbulence have a lower luminosity. 



7 METAL LINES 

Turbulent gas motions can lead to substantial Doppler broadening 
of X-ray metal lines, in excess of their intrinsic line widths. Given 
the exquisite spectral resolution of upcoming observational X-ray 
mission, this could be used to directly measure the degree of ICM 
turbulence (Sunvaev et aL .2003) by measuring, e.g., the shape of 
the 6.702 keV iron emission line. 

One potential difficulty in this is that multiple large-scale bulk 
motions of substructure moving inside the galaxy cluster along the 
line of sight might dominate the signal. To get a handle on this, 
we estimate the line-of-sight emission of the 6.702 keV iron line 
within columns through the simulated clusters, where the column 
size was chosen to be 300 /i~^kpc on a side, which at the distance 
of the Coma cluster corresponds roughly to one arcmin, the formal 
resolution of ASTRO-E2. For simplicity, we assign every gas par- 
ticle a constant iron abundance and an emissivity proportional to 
nl X /(Te) X AV, where Ue is the electron density, and AV (x 
is the volume represented by the particle. As a further approxima- 




Figure 13. Radial velocity dispersion profile for dark matter (black) and 
gas (blue) particles. The thick lines represent the average over the 4 massive 
clusters, whereas the thin lines give the average over the 5 low mass sys- 
tems. The dashed lines are drawn from the original viscosity simulations, 
the solid lines from the low-viscosity simulations. 
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Figure 14. Radial gas density profile. The thick lines represent the average 
over the 4 massive clusters, whereas the thin lines give the average over the 
5 low mass systems. The dashed lines are drawn from the original viscosity 
simulations, the solid lines from the low-viscosity simulations. 
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Figure 15. Mass -weighted gas temperature profile. The thick lines repre- 
sent the average over the 4 massive clusters, whereas the thin lines give the 
average over the 5 low mass systems. The dashed lines are drawn from the 
original viscosity simulations, the solid lines from the low-viscosity simu- 
lations. 
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Figure 16. Radial entropy profiles of the ICM gas. Thin lines are individual 
profiles for the 9 clusters, thick lines are averages. The dashed lines are 
drawn from the original viscosity simulations, the solid lines from the low- 
viscosity simulations. 



tion we set the electron abundance equal to unity. We also neglect 
thermal broadening and other close lines (like the 6.685 keV iron 
line), given that the 6.702 keV iron line is clearly the strongest. 

In Figure fTsl we show the resulting distributions for several 
lines of sight, here distributed on a grid with —500, —250, 0, 250 
and 500 /i~^kpc impact parameter in x-direction, and —250, and 
250, h~^kpc impact parameter in the ^/-direction, respectively. The 
different lines in each panel correspond to simulations with the 
signal-velocity based viscosity (dashed line) and with the time- 
depended low-viscosity scheme (solid lines). Both results have 
been normalised to the total cluster luminosity, such that the in- 
tegral under the curves corresponds to the fraction of the total lu- 
minosity. 

We note that this measurement is very sensitive to small tim- 
ing differences between different simulations, and therefore a com- 
parison of the same cluster run with different viscosity should be 
carried out in a statistical fashion, even if some features look very 
similar in both runs. In gene ral we confirm previous findings (e.g. 
Ilnogamov & Sunvaevll2003h that large bulk motions can lead to 
spectral features which are several times wider than expected based 
on thermal broadening alone. Additional complexity is added by 
beam smearing effects, thermal broadening, and by the local turbu- 
lence in the ICM gas, such that an overall very complex line shape 
results. In our simulations with the low-viscosity scheme, where 
we have found an increased level of fluid turbulence, the final line 
shapes are indeed more washed out. However, the complexity of the 
final line shapes suggests that it will be very difficult to accurately 
measure the level of fluid turbulence with high resolution spectra of 
X-ray emission lines, primarily because of the confusing influence 
of large-scale bulk motions within galaxy clusters. 
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Figure 17. Distribution of the distance of particles to their nearest halo at high redshift. The particles selected here end up within 5% of i?vir at z = 0. The 
dashed lines are for the original viscosity scheme, while the solid lines mark the result for the low-viscosity simulations. 



8 APPLICATION TO RADIO HALOS 

One promising possibility to explain the extended radio emission 
on Mpc- scales observed in a growing number of galaxy clusters is 
to attribute it to e lectron acceleration by cluster turbulence (e.g., 
ISchlickeiser et alJf lQSV: Brunetti et al. 2001). Having high resolu- 
tion cluster simulations at hand, which thanks to the new viscosity 
scheme are able to develop significant turbulence within the ICM, 
it is of interest to explore this possibility here. Obviously, due to the 
uncertainties in the complex physical processes of dissipation of the 
turbulent energy - which up to this point can not be explicitly mod- 
elled in the simulations - our analysis is limited to a check whether 
or not turbulent reacceleration can plausibly reproduce some of the 
main properties of radio halos. In this scenario, the efficiency of 
electron acceleration depends on the energy density of magneto- 
hydrodynamic waves (Alfven waves. Fast Mode waves, . . . ), on 
their energy spectrum, and on the physical conditions in the ICM 
(i.e., density and temperature of the thermal plasma, strength of the 
magnetic field in the ICM, number density and spectrum of cosmic 
rays in the ICM). A number of approaches for studying the accel- 
eration of relativistic electrons in the ICM have been succe ssfully 
devel oped by focusing on the case of Alfven waves ( Oh no et alJ 

f?;'Brunetti et al."2004') and, more recently, on Fast Mode- waves 
ssano & Brunetti 2005). 
It should be stressed, however, that analytical and/or semi- 
analytical computations are limited to very simple assumptions for 
the generation of turbulence in the ICM. Full numerical simulations 
represent an ideal complementary tool for a more general analysis, 
where the injection of turbulence into the cluster volume by hierar- 
chical merging processes can be studied realistically. Low numeri- 
cal viscosity and high spatial resolution are however a prerequisite 
for reliable estimates of turbulence. As we have seen earlier, previ- 
ous SPH simulations based on original viscosity parameterisations 
have suppressed random gas motions quite strongly, but the low- 
viscosity scheme explored here does substantially better in this re- 
spect. 

In this Section, we carry out a first exploratory analysis of 
the efficiency of electron acceleration derived in the low-viscosity 
scheme, and we compare it to results obtained with a original SPH 
formulation. For definiteness, we assume that a fraction rjt of the 
estimated energy content of the turbulent velocity fields in the clus- 
ter volume, measured by the local velocity dispersion (equation 
[T6l section |5}, is in the form of Fast Mode waves. We focus on 



these modes since relativistic electrons are mainly accelerated by 
coupling with large scale modes (e.g., ^ kpc, k being the 
wave number) whose energy density, under the above assumption, 
can hopefully be constrained with the numerical simulations in a 
reliable fashion. In addition, the damping and time evolution of 
Fast Modes basically depend only on the properties of the thermal 
plasma and are essentia lly not sensitive t o the presence of cosmic 
ray protons in the ICM (Cassano & Brunetti 2005). 

Relativistic particles couple with Fast Modes via magnetic 
Landau damping. The nec essary condition for Landau damping 
d Melrosel [1963: lEilekl ITqTq ) is cj - A;||^;|| 0, where uj is the 
frequency of the wave, /cy is the wavenumber projected along the 
magnetic field, and v\\ = t'/i is the projected electron velocity. Note 
that in this case - in contrary to the Alfvenic case - particles may 
also interact with large scale modes. In the collisionless regime, it 
can be shown that the resulting acceleration rate in an isotropic 
plas ma (modes' propagation a nd particle momenta) is given by 
(e.g.,|Cassano & Brunett|2005|) 



(23) 



where vm is the magneto-sonic velocity, and is the energy 
spectrum of the magnetic field fl uctuations re.g.. lBarnes & Scargld 
|1273|; [Cassano & Brunetti||200j. We estimate the rate of injection 
of Fast Modes, ^, assuming that a fraction, rft, of fluid turbu- 
lence is associated with these modes. We parameterise the injection 
rate assuming that turbulence is injected (and also dissipated) in 
galaxy clu sters within a time of the order of a cluster-cros sing time, 
Tcross fsee Cassano & Brunetti[l2005l: [Subramanian et"al]r2005L for 
a more detailed discussion). One then has: 



/ 



jFM 1 7 



' rit- 



Et 



1 2 -1 



(24) 



Here, Tinject is the time over which a merging substructure 
is able to inject turbulence in a given volume element in the main 
cluster. This can be estimated as the size of the subhalo divided by 
its infalling velocity. As the size of a halo is only a weak function of 
its mass, we approximate Tinject with a generic value of Tinject = 
0.5 Gyr. This is only a very crude estimate and more generally one 
should think of an effective efficiency parameter 77?^ — T\tl Tinject 
which we set to 0.1/ (0.5 Gyr) as argued before. Note also that for 
estimating gI we used a 64^ TSC-grid, which is a conservative 
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Figure 18. Distribution of the Doppler-shifted emission of the iron 6.702 keV line for 15 lines of sight through the cluster g72. Every panel corresponds to a 
column of side length 300 h~^kpc through the virial region of the cluster. This roughly corresponds to one arcmin resolution (comparable to the ASTRO-E2 
specifications) at the distance of the Coma cluster. The columns from left to right correspond to —500, —250, 0, 250, and 500 h~^kpc impact parameter 
in the x-direction, and the rows correspond to —250, 0, and 250, /i~^kpc impact parameter in the y-direction. The dashed lines give results for the original 
viscosity run, while the solid line is for the low-viscosity run. The thick bar in the center panel marks the expected energy resolution of 12 eV of the XRS 
instrument on-board ASTRO-E2. 



estimate, as shown in Figure |4| and therefore equation l24l should 
still reflect a lower limit. 

Following Cassano & Brunetti (2005), the spectrum of the 
magnetic fluctuations associated with Fast Modes is computed un- 
der stationary conditions taking into account the damping rate of 
these modes with thermal electrons, Tk = Tok. One then has 



^2 



1 Jf^ 



8n Pgas Tok 



(25) 



Thus the integral in Eqn. I23> at each position of the grid can be 
readily estimated as 



/ 



kWk dk 



B. 
Stt 



i 1 f 

Fo-Pgas J 



tFM 1 7 

iz, dk 



eff ^^(x) Pgas(x)cr,^i(x) 1 



IGtt 



Pgas (x) Fo 



(26) 
(27) 



where Fo depends on the temperature of the ICM (Cassano & 
Brunetti 2005) \ 

In this Section we are primarily interested in determining the 
maximum energy of accelerated electrons, given the adopted en- 
ergy density for Fast Modes. Under typical conditions in the ICM, 
the maximum energy of electrons is reached at energies where ra- 
diative losses balance the effect of the acceleration. The radiative 
synchrotron a nd inver se Compton losses with CMB photons are 
given by fe.g.JSarazinfclQQa) 



(I) 



-4.8 X 10" V 



2/3 



+ il + z) 



rUe c' 



(28) 



^ Note that under these assumptions the efficiency of the particle accelera- 
tion does not depend on the spectrum of the waves 
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where B^g is the magnetic field strength in jiG, and is the pitch 
angle of the emitting electrons. If an efficient isotropisation of elec- 
tron momenta can be assumed, it is possible to average over all 
possible pitch angles, so that ^sin^ 6^ — 2/3. 

In Figure fTol we plot the maximum energy of the fast electrons 
obtained from Eqs. J25j and J28j along one line-of- sight through 
the cluster atmosphere. The two different lines are for the same 
cluster, simulated with our two main schemes for the artificial vis- 
cosity. The two vertical lines in Figure|5]are indicating the position 
of these cuts. When the new low- viscosity scheme is used, enough 
turbulence is resolved to maintain high energy electrons (and thus 
synchrotron radio emission) almost everywhere out to a distance of 
1 Mpc from the cluster centre, whereas in the original formulation 
of SPH, turbulence is much more strongly suppressed, so that the 
maximum energy of the accelerated electrons remains a factor of 
about three below that in the low-viscosity case. 

The results reported in Figure[T9lare obtained assuming a ref- 
erence value of rjif^ — r]t /rinject — 0.1/(0.5 Gyr). The averaged 
volume weighted magnetic field strength in the considered cluster 
region is fixed at 0.5/j.G and a simple scaling from magnetic flux- 
freezing, B adopted in the calculations, resulting in 
a central magnetic field strength of Bo = 5.0/j.G. It is worth not- 
ing that the maximum energy of the accelerated electrons, 7max, 
scales with the energy density of the turbulence (and with the frac- 
tion of the energy of this turbulence channelled into Fast Modes rit, 
EQ. I25t . However the synchrotron frequency emitted by these elec- 
trons scales with the square of the turbulent energy (7max)- Interest- 
ingly, with the parameter values adopted in Figure[l9] a maximum 
synchrotron frequency of order 10^(?7t/0.1)^ MHz is obtained in 
a Mpc-sized cluster region, which implies diffuse synchrotron ra- 
diation up to GHz frequencies if a slightly larger value of rjt is 
adopted^ . On the other hand, we note that essentially no significant 
radio emission would be predicted if we had used the simulations 
with the original SPH viscosity scheme. 

In real galaxy clusters, the level of turbulence which can form 
will also depend on the amount of physical viscosity present in the 
ICM gas (i.e. on its Reynolds number), which in turn depends on 
the magnetic field topology and gas temperature. It will presumably 
still take a while before the simulations achieve sufficient resolution 
that the numerical viscosity is lower than this physical viscosity. In 
addition, the details of the conversion process of large-scale ve- 
locity fields into MHD modes is still poorly understood and well 
beyond the capabilities of presently available cosmological simula- 
tions. However, our results here show that a suitable modification 
of the artificial viscosity parameterisation within SPH can be of sig- 
nificant help in this respect, and it allows a first investigation of the 
role of turbulence for feeding non-thermal phenomena in galaxy 
clusters. 



9 CONCLUSIONS 

We implemented a new parameterisation of the artificial viscos- 
ity of SPH in the parallel cosm ological simulation code GA DGET- 
2. Following a suggestion by iMorr is & Monaghai] ll991 ). this 
method amounts to an individual, time-dependent strength of the 
viscosity for each particle which increases in the vicinity of shocks 
and decays after passing through a shock. As a result, SPH should 

^ Note that re.g.. lCassano & Brunettil2005l) required r]t ~ 0.2 - 0.3 in 
order to reproduce the observed statistics of radio halos. 
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Figure 19. One-dimensional profile of the maximum energy of the electrons 
accelerated via the turbulent-magneto-sonic model, along the same vertical 
lines drawn in Figure|5] Dashed lines are for the original viscosity run, while 
solid lines are for the low-viscosity scheme. Here, a conservative 64^ grid 
is used in the TSC smoothing. 

show much smaller numerical viscosity in regions away from 
strong shocks than original formulations. We applied this low- 
viscosity formulation of SPH to a number of test problems and to 
cosmological simulations of galaxy cluster formations, and com- 
pared the results to those obtained with the original SPH formula- 
tion. Our main results can be summarised as follows: 

• The low-viscosity variant of SPH is able to capture strong 
shocks just as well as the original formulation, and in some cases 
we even obtained improved results due to a reduced broadening 
of shock fronts. In spherical accretion shocks, we also obtained 
slightly better results due to a reduction of pre- shock entropy gen- 
eration. 

• Using the low-viscosity scheme, simulated galaxy clusters 
developed significant levels of turbulent gas motions, driven by 
merger events and infall of substructure. We find that the kinetic 
energy associated with turbulent gas motion within the inner ^ 
1 Mpc of a 10^^ h~^MQ galaxy cluster can be up to 30% of the 
thermal energy content. This value can be still larger and reach up 
to 50% in the very central part of massive clusters. In clusters with 
smaller masses 10^^ /i^^M©) we find a smaller turbulent en- 
ergy content, reaching only 5% within the central Mpc. Within a 
comparable fraction of the virial radius, the corresponding fraction 
is however still of order 10%. These values are much larger than 
what is found when the original SPH viscosity is employed, which 
strongly suppresses turbulent gas motions. 

• The presence of such an amount of turbulence has an imprint 
on global properties of galaxy clusters, most notably reducing the 
bolometric X-ray luminosity in non radiative simulations by a fac- 
tor of ^ 2. However, the global, mass- weighted temperature does 
not change. 

• The temperature profiles of galaxy clusters are only mildly 
changed by the presence of turbulence, but we observe a strong de- 
crease of density within the central region of galaxy clusters, where 
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the turbulence is providing a significant contribution to the total 
pressure. Also the radial entropy profiles show a significant flatten- 
ing towards the cluster centre. This makes them more similar to the 
observed profiles based on X-ray observations. Note however that 
radiative cooling - which was not included in our simulations - can 
also modify the profiles substantially. We find that the higher en- 
tropy in the centre found in the low viscosity simulations is largely 
a result of the more efficient transport and mixing of low-entropy 
in infalling material into the core of the cluster. We note that the 
elevated entropy levels found in our low-viscosity runs are more 
similar to the results found with Eulerian hydrodynamic codes than 
the original SPH ones. 

• Turbulence in galaxy clusters broadens the shape of metal 
lines observable with high-resolution X-ray spectrographs like 
XRT on board of ASTRO-E2. Depending on the strength of the tur- 
bulence and the dynamical state of the cluster, prominent features 
due to large-scale bulk motions may however get washed out and 
blended into a very complex line structure. In general it will there- 
fore be difficult to isolate the signature of the turbulent broadening 
and to differentiate it unambiguously from the more prominent fea- 
tures of large scale bulk motions. 

• Applying a model for accelerating relativistic electrons by 
ICM turbulence we find that galaxy clusters simulated with reduced 
viscosity scheme may develop sufficient turbulence to account for 
the radio emission that is observed in many galaxy clusters, pro- 
vided that a non-negligible fraction of the turbulent energy in the 
real ICM is associated with Fast Modes. 

In summary, our results suggest that ICM turbulence might be 
an important ingredient in the physics of galaxy clusters. If present 
at the levels inferred from our low-viscosity simulations, it has a 
significant effect on the radial structure and on the scaling relations 
of galaxy clusters. We also note that the inferred reduction of the 
X-ray luminosity has a direct influence on the strength of radiative 
cooling flows. The more efficient mixing processes would also help 
to understand the nearly homogeneous metal content observed for 
large parts of the cluster interior. Finally, cluster turbulence may 
also play an important role for the dynamics of non-thermal pro- 
cesses in galaxy clusters. 

Although we observe a rather high level of turbulence in the 
very centre of our simulated galaxy clusters when we use the low- 
viscosity scheme, it is likely that we are still missing turbulence 
due to the remaining numerical viscosity of our hydrodynamical 
scheme, and due to the resolution limitations, particularly in low 
density regions, of our simulations. This problem should in prin- 
ciple become less and less severe as the resolution of the simula- 
tions is increased in future calculations. However, given that there 
is a some physical viscosity in real galaxy clusters which limits the 
Reynolds number of the ICM, it cannot be the goal to model the 
ICM as a pure ideal gas. Instead, future work should concentrate on 
accurately characterising this physical viscosity of the ICM, which 
could then be directly incorporated into the simulations by means 
of the Navier-Stokes equations. Our results suggest that the low- 
viscosity formulation of SPH should be of significant help in re- 
ducing the numerical viscosity of SPH simulation below the level 
of this physical viscosity, and the present generation of simulations 
may already be close to this regime. 
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